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It  is  the  purpose  of  this  dissertation  to  investi- 
gate  the  dynamic  stability  of  elastic  shells  of  revolu- 
tion.  Two  specific  areas  of  this  broad  field  are  treated 
in  detail.   First,  this  analytical  study  generates  pre- 
viously unavailable  interaction  relations  for  combined 
dynamic  loadings  as  they  interact  to  cause  passage  from  a 
dynamically  stable  state  to  other  states.   Secondly,  the 
concept  of  linearized  dynamic  stability  is  extended  to 
include  geometrically  nonlinear  effects.   One  role  of  these 
effects  is  to  allow  the  possibility  of  autoparametri c  ex- 
citation of  preferred  flexural  modes  by  the  driven  modes. 
The  question  of  whether  such  excitation  can  occur  during 
the  primary  dynamic  instability  motion  sufficiently  to  af- 
fect the  magnitude  of  critical  dynamic  loads  is  studied. 


A  modified  version  of  the  subdomain  method  in  combination 
with  circumferential  modal  analysis  is  developed  for  the 
solution  method.   A  computer  program  is  constructed  to 
obtain  numerical  results  and  the  dynamic  stability  charac- 
teristics of  a  conical  frustum  are  studied  under  a  variety 
of  combined  dynamic  loadings.   Interaction  curves  for 
static  stability,  linearized  dynamic  stability,  and  non- 
linear dynamic  stability  are  generated. 

For  the  particular  loading  conditions  studied,  the 
results  indicate  that  dynamic  instability  will  occur  at 
loads  below  critical  static  loads.   This  reduction  in  crit- 
ical dynamic  loads  was  shown  to  be  the  result  of  both  the 
dynamic  load  factor  effect  and  of  autoparametri c  excitation. 
The  interaction  curves  which  are  generated  illustrate  these 
effects  quantitatively  under  conditions  of  combined  dynamic 
loading.   A  criterion  for  dynamic  buckling  is  established 
based  on  meridional  mode  shape  changes.   This  ability  to 
detect  sudden  jumps  in  the  meridional  profile  helps  to  verify 
instability  detected  by  divergence  on  displacement  time 
history  curves  and  provides  additional  information  about 
the  poststable  state. 


XT 


CHAPTER  1 


INTRODUCTION 


A  dynamic  stability  problem  can  be  defined  as  a 
problem  analyzed  by  Newton's  equations  of  motion  or  any 
equivalent  method  [1].   Dynamic  stability  problems  for 
continuous  systems  are  nearly  always  governed  by  nonlinear 
partial  differential  equations  [2].   The  solution  to  a 
problem  of  dynamic  stability  in  thin-walled  structures  may 
be  stated  as  the  determination  of  what  load-time  combina- 
tions cause  the  displacements  of  points  in  the  material 
body  to  increase  sufficiently  to  either  interfere  with 
operational  specifications  or  cause  a  breakdown  of  the 
structure.   Methods  for  obtaining  such  solutions  may  be 
categorized  into  three  fundamental  areas.   The  most  obvi- 
ous and  usually  most  difficult  method  is  the  direct  inte- 
gration of  the  governing  equations  for  sufficiently  long 
periods  to  allow  direct  observation  of  the  ensuing  motions. 
A  second  technique  is  the  investigation  of  the  stability 
of  equilibrium  points.   The  question  here  is  whether  a 
slight  perturbation  of  a  dynamical  system  from  an  equilib- 
rium state  will  produce  motion  confined  to  the  neighborhood 


of  the  equilibrium  point  or  a  motion  tending  to  leave 
that  neighborhood  [3].   This  method  is  usually  based  on 
either  perturbation  or  characteristic  equation  techniques. 
A  third  category  of  dynamic  stability  methods  is  on  a 
higher  level  of  abstraction  than  the  previous  two  and 
relies  on  the  existence  and  properties  of  variously  de- 
fined functionals.   Energy  functionals  [4]  and  Liapunov 
functionals  [5]  are  two  examples. 

Definitions  of  dynamic  stability  may  be  categorized 
by  methods  of  approach  to  such  problems.   In  situations 
where  the  direct  integration  approach  is  employed,  dynamic 
instability  may  be  defined  to  exist  when  some  character- 
istic parameter,  such  as  deflection,  stress,  or  strain, 
becomes  unbounded  in  one  or  more  parts  of  the  body.   When 
the  method  of  small  oscillations  about  equilibrium  points 
is  used,  instability  may  be  defined  as  unbounded  growth 
of  perturbations  with  time  or  by  the  existance  of  one  or 
more  positive  real  parts  of  roots  of  a  characteristic 
equation.   In  problems  where  the  solution  method  is  of 
the  third  above  mentioned  category  a  sufficient  condition 
for  stability  is  defined  according  to  the  definite  form 
of  the  functional  employed.   In  the  case  of  energy  func- 
tionals, a  boundary  between  stable  and  unstable  regions 
is  defined  by  the  vanishing  of  the  second  variation  of 
the  energy  functional.   When  Liapunov  functionals  are 


employed,  stability  is  defined  when  the  functional  is 
positive  definite  and  when  its  derivative  is  negative 
def i  ni  te  . 

The  particular  problem  of  interest  in  this  study 
is  the  dynamic  buckling  of  shell  structures,  which  rep- 
resents a  specific  application  of  the  more  general  theory 
of  dynamic  stability  principles.   Just  as  the  definitions 
of  dynamic  stability  are  related  to  the  solution  methods, 
the  definitions  of  dynamic  buckling  of  shells  are  related 
to  solution  methods  as  well  as  the  specifics  of  the  phys- 
ical problem  such  as  geometry,  material  properties,  or 
loading  rates. 

Solution  methods,  falling  within  the  first  of  the 
three   aforementioned  categories,  are  direct  temporal  in- 
tegration of  equations  of  motion.   They  may  be  classified 
into  five  procedural  divisions.   These  are  finite  differ- 
ence, finite  element,  Hamilton's  principle,  energy  methods, 
and  methods  of  weighted  residuals.   Figure  1  illustrates 
the  procedural  flow  and  interrelations  between  these  meth- 
ods for  the  study  of  dynamic  stability  of  elastic  shells 
and  is  based  on  the  most  prominant  articles  dealing  with 
this  subject. 

First,  let  us  consider  recent  significant  contri- 
butions falling  within  the  finite  difference  method  divi- 
sion.  Witmer  et  al.  [6,  7,  8]  have  developed  during  the 
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Figure  1.   Development  Paths  in  Constructing  Solution 
Methods  for  Problems  of  Dynamic  Stability 
of  Shells 


past  decade  a  shell  dynamics  computer  program  known  as 
PETROS  which  utilizes  spatial  finite  difference  methods 
in  both  the  circumferential  and  meridional  directions. 
Utilizing  PETROS  the  authors  have  studied  a  number  of 
explosive  forming  problems  and  have  shown  remarkable 
agreement  with  experimental  results.   However,  consider- 
ing that  after  a  ten-year  development  period,  computer 
run  times  are  still  measured  in  hours  for  moderately  com- 
plex problems,  the  practical  usefulness  of  PETROS  in  study- 
ing shell  dynamic  stability  is  somewhat  limited.   Using 
shallow  shell  theory,  Longhitano  [9]  developed  numerical 
solutions  for  the  nonlinear  dynamics  of  hemispherical  shells 
also  utilizing  finite  difference  solutions  in  two  directions 
He  showed  that  suddenly  applied  symmetic  pressures  may 
couple  the  breathing  mode  into  more  preferred  axisymmetric 
or  asymmetric  vibratory  modes  causing  instability  at  approx- 
imately one  half  static  buckling  pressures.   Again,  however, 
the  two-dimensional  discrete  type  of  spatial  treatment  re- 
sults in  unusually  long  computer  run  times  thus  limiting 
the  number  of  cases  treated. 

Another  class  of  finite  difference  solutions  to 
shell  dynamic  stability  problems  have  appeared  in  the  lit- 
erature which  require  only  approximately  three  fourths  of 
the  computer  run  times  as  the  previously  cited   two-dimen- 
sional discrete  methods  and  appear  to  give  comparable 


results.   Such  solutions  based  on  elimination  of  the 
circumferential  coordinate  via  modal  analysis  prior  to 
finite  differencing  in  the  meridional  direction  have 
been  presented  by  Cromer  [10],  Ball  [11,  12],  and  Kim 
[13].   The  results  of  Cromer  and  Ball  also  indicate  that 
suddenly  applied  axisymmetric  pressures  may  couple  breath- 
ing mode  response  into  preferred  vibratory  modes  to  cause 
dynamic  instabilities  at  pressure  loadings  lower  than 
static  critical  loads.   These  studies  indicate  that  cou- 
pling between  parametrically  excited  flexural  modes  is 
very   weak  and  often  may  be  neglected  without  compromising 
the  resul ts . 

Stricklin  [14]  and  Wu  [15]  have  made  recent  and 
significant  contributions  to  the  nonlinear  dynamic  analy- 
sis of  shells  by  employing  the  finite  element  method.   Re- 
sults published  by  both  of  these  investigators  indicate 
that  accuracy  comparable  to  methods  emanating  from  the 
equations  of  motion  is  obtainable  if  sufficiently  complex 
elements  are  used  in  large  quantity.   However,  due  to 
longer  computing  times  required  for  finite  element  tech- 
niques, their  major  advantage  clearly  lies  in  the  treat- 
ment of  complex  structures  where  equations  of  motion  are 
not  reasonably  derivable.   A  seldom  mentioned,  yet  almost 
universal,  appealing  characteristic  of  finite  element 
techniques  warrants  a  brief  mention.   This  is  that  the 


unknown  dependent  variables  in  a  finite  element  solution 
have  a  direct  physical  interpretation  at  all  times  as 
the  solution  proceeds.   Thus,  to  some  extent,  the  analyst 
never  loses  physical  contact  with  his  problem.   This  is 
in  contrast  to  most  other  analysis  techniques  applicable 
to  shell  dynamics  where  the  solution  is  carried  forward 
in  some  transformed  mathematical  space. 

Klosner  [16,  17]  in  a  series  of  studies  covering 
a  five-year  period  studied  the  effect  of  suddenly  applied 
axial  loads  on  cylindrical  shells.   His  solutions  were  ob- 
tained by  elimination  of  the  circumferential  coordinate 
from  nonlinear  shallow  shell  equations  and  applying  the 
Ritz  method  or  Hamilton's  principle  to  the  relations  re- 
sulting after  introducing  assumed  solutions.   His  results 
indicate  that  suddenly  applied  axial  loads  do  not  exhibit 
the  dynamic  load  factor  effect  found  by  others  [9,  10,  11, 
12]  when  suddenly  applied  normal  pressures  are  considered. 
These  investigations  reveal  that  a  valuable  indication  of 
the  onset  of  dynamic  instability  can  be  found  in  the  re- 
versal in  the  trend  of  the  time  to  first  maximum  deflec- 
tion of  the  critical  mode. 

Another  subcategory  of  solution  methods  in  studying 
dynamic  shell  stability  by  direct  temporal  integration  is 
known  as  the  energy  method.   This  technique  circumvents 
the  derivation  of  the  partial  differential  equations  of 


motion.   The  method  involves  deriving  internal  energy 
expressions  in  conjunction  with  appropriate  kinematic 
relations,  substituting  assumed  solutions  into  these 
relations,  and  then  applying  Lagrange's  equations  to 
the  result.   This  procedure  leads  to  either  a  coupled 
system  of  nonlinear  ordinary  differential  equations  if 
the  spatial  integrals  are  evaluated  analytically  or  a 
coupled  system  of  nonlinear  i ntegro-di f f erenti al  equa- 
tions if  the  spatial  integrals  are  to  be  evaluated 
numerically.   In  either  case,  the  differentials  are 
temporal  and  must  be  solved  numerically.   Goodier  [18], 
Mclvor  [19],  and  Lindberg  [20]  have  presented  solutions 
to  a  particular  type  of  dynamic  instability  of  cylin- 
drical shells  utilizing  this  solution  procedure.   Hubka 
[21]  has  presented  alternate  methods  for  numerical  inte- 
gration of  the  final  equation  sets  for  this  problem. 
This  problem  is  concerned  with  the  circumferential  modal 
coupling  arising  from  the  nonlinear  membrane  stiffness 
terms  in  the  equation  sets  obtained  by  this  energy  method 
technique.   The  main  results  of  these  studies  are  that  non- 
linear coupling  between  circumferential  modes  may  under 
conditions  of  minor  excitation  or  initial  imperfections 
in  the  preferred  vibratory  modes  cause  energy  extraction 
from  the  driven  breathing  mode  into  one  or  more  preferred 
modes  such  that  the  preferred  mode  becomes  dynamically 
unstabl e . 


utilizing  this  same  energy  method  derivation, 
Mente  [22]  has  studied  the  dynamic  nonlinear  response 
of  cylindrical  shells  subjected  to  asymmetric  pressure 
loading.   His  computerized  solution,  known  as  DEPICS, 
is  able  to  handle  a  variety  of  coupled  modes  during  the 
temporal  integration.   The  nonlinear  stiffness  is  handled 
by  solving  a  linear  eigenvalue-eigenvector  problem  at 
each  time  step  then  adjusting  the  succeeding  stiffness 
value  based  on  the  nonlinear  terms.   Unfortunately  this 
rather  extensive  computer  program  has  to  date  been  unable 
to  satisfactorily  define  any  regions  of  instability. 

Significant  Russian  contributions  to  the  study  of 
dynamic  stability  of  shells  have  been  made.   Among  these 
are  works  by  Shumik  [23,  24]  who  also  employed  solutions 
by  direct  temporal  integration  of  expressions  derived  by 
the  energy  method.   These  works  deal  with  suddenly  applied 
loads  using  a  linearized  uncoupled  dynamic  stability  theory 
His  results  indicate  that  cones  with  large  taper  buckle 
dynamically  at  circumferential  wave  numbers  higher  than 
the  critical  wave  number  for  static  loading. 

Shiau  [25]  utilizing  an  energy  solution  studied 
the  effects  of  imperfections  on  the  dynamic  stability  of 
conical  shells  subject  to  axial  impact.   His  results 
agree  with  Klosner's  [16,  17]  in  so  much  that  a  reversal 
in  the  trend  of  the  time  to  first  maximum  deflection  of 
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the  critical  mode  may  be  viewed  as  an  indication  of  the 
onset  of  dynamic  buckling.   The  results,  however,  are 
so  critically  dependent  upon  the  choice  of  the  single  term 
assumed  buckling  mode  that  the  procedure  appears  to  be 
of  limited  general  value. 

A  final  division  of  solution  methods  in  the  cate- 
gory of  direct  temporal  integration  is  the  method  of 
weighted  residuals.   The  details  of  the  procedure  will 
be  covered  in  a  later  chapter.   Nash  [26]  utilizing 
shallow  shell  theory  in  conjunction  with  the  Galerkin 
method  of  weighted  residuals  investigated  the  response 
of  thin  conical  shells  to  dynamically  applied  axial 
forces.   He  found  that  a  rapid  increase  in  end  shorten- 
ing may  be  used  to  denote  the  onset  of  dynamic  instability. 

Fisher  [27]  employed  a  similar  solution  method  to 
study  the  dynamic  buckling  of  reinforced  circular  cylin- 
ders subjected  to  suddenly  applied  axial  compression  load- 
ings.  He  also  found  that  a  reversal  of  the  trend  of  the 
time  to  maximum  deflection  corresponds  to  onset  of  dynamic 
buckling.   Lakshmikantham  [28]  using  shallow  shell  theory 
in  conjunction  with  smeared  analysis  and  the  Galerken  tech- 
nique investigated  the  same  problem  as  Fisher  [27].   His 
solutions  indicated  that  a  jump  in  radial  deflection  was 
the  correct  criteria  for  indicating  dynamic  instability. 
An  earlier  work  by  Dietz  [29]  studying  the  same  problem 
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by  shallow  shell  theory  and  the  Galerkin  method  indi- 
cated that  loss  of  stability  is  best  defined  as  a  sudden 
drop  in  the  axial  load  required  to  maintain  a  constant 
end  shortening  rate. 

By  taking  the  results  of  the  preceding  literature 
survey  in  context,  several  comments  may  be  made  concerning 
the  current  state-of-the  art  in  analytical  studies  of 
dynamic  shell  stability.   First,  it  appears  that  nonlinear 
response  problems  which  retain  full  circumferential  modal 
coupling  yield  the  most  reliable  indicators  of  the  onset 
of  dynamic  buckling.   Also  the  results  appear  to  be  highly 
sensitive  to  theoretical  approximations  and  to  the  quality 
of  the  numerical  solution  procedures.   Secondly,  the  study 
of  passage  through  bifurcation  points  is  well  developed 
only  for  linearized  dynamic  stability  techniques.   The 
few  articles  proposing  to  deal  with  nonlinear  dynamic 
stability  unfortunately  discarded  the  quadratic  nonline- 
arities,  which  are  the  most  significant  ones,  retaining  in- 
stead only  the  cubic  nonl i neari ti es  .   This  was  done  in  order 
to  circumvent  having  to  deal  with  inter-modal  coupling.  There- 
fore, the  results  of  these  studies  appear  inconsistent. 
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It  is  the  purpose  of  this  dissertation  to  invest- 
igate the  dynamic  stability  of  elastic  shells  of  revolu- 
tion.  Two  specific  areas  of  this  broad  field  are  treated 
in  detail.   First,  this  analytical  study  generates  previ- 
ously unavailable  interaction  relations  for  combined 
dynamic  loadings  as  they  interact  to  cause  passage  from 
a  dynamically  stable  state  to  other  states.   Secondly,  the 
concept  of  linearized  dynamic  stability  is  extended  to  in- 
clude geometrically  nonlinear  effects.   One  role  of  these 
effects  is  to  allow  the  possibility  of  autoparametri c 
excitation  of  preferred  flexural  modes  by  the  driven  modes. 
The  question  of  whether  such  excitation  can  occur  during 
the  primary  dynamic  instability  motion  sufficiently  to  af- 
fect the  magnitude  of  critical  dynamic  loads  is  studied. 

The  theory  derived  herein  is  suitable  for  the  study 
of  the  dynamic  stability  of  general  shells  of  revolution. 
In  order  to  obtain  numerical  results,  however,  a  specific 
shape  must  be  specified  and  the  one  used  for  numerical 
analysis  in  this  study  is  a  truncated  circular  conical 
shell  structure.   Conical  frustum  shells  are  in  frequent 
use  as  structural  elements.   Moreover,  a  conical  shell 
reduces  to  a  cylindircal  shell  when  the  semivertex  angle 
of  the  cone  becomes  zero.   Similarly,  it  reduces  to  a  flat 
circular  plate  when  the  semivertex  angle  of  the  cone  ap- 
proaches a  right  angle.   Thus,  the  analysis  and  computer 
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program  developed  for  cones  can  be  generalized  to  handle 
these  two  related  problems. 

A  modified  version  of  the  subdomain  method  in 
combination  with  circumferential  modal  analysis  is  de- 
veloped for  the  solution  method.   A  computer  program  is 
constructed  to  obtain  numerical  results  and  the  dynamic 
stability  characteristics  of  a  conical  frustum  are  stud- 
ied under  a  variety  of  combined  dynamic  loadings.   Inter- 
action curves  for  static  stability,  linearized  dynamic 
stability,  and  nonlinear  dynamic  stability  are  generated, 
A  dynamic  buckling  criterion  based  on  meridional  mode 
shape  changes  is  established  and  the  quantitative  impor- 
tance of  nonlinear  coupling  effects  is  illustrated. 


CHAPTER  2 
THEORETICAL  DEVELOPMENT  OF  NONLINEAR  SHELL  EQUATIONS 

2  . 1   Shell  Coordinates 

A  set  of  normal  curvilinear  coordinates  (a,  6,  z) 
is  chosen  in  the  shell  space  such  that  a  and  6  lie  on  the 
undeformed  middle  surface,  z  is  perpendicular  to  this 
middle  surface,  and  z-0  lies  at  the  middle  surface.   The 
equation  of  the  undeformed  middle  surface  is  given  in 
terms  of  two  independent  parameters  a  and  B  by  the  radius 
vector 


r  =  r(a,3) . 


2.1) 


To  describe  the  location  of  an  arbitrary  point  in  the 
space  occupied  by  a  shell,  the  position  vector  is  defined 
as 


R(a,3,z)  =  r(a,6)  +  zi  (a, 6) 


(2.2) 


where  z  measures  the  distance  of  the  point  from  the  corres 
ponding  point  on  the  middle  surface  along  i  .   The  unit 
vector  i^  is  normal  to  the  surface. 
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2.2   Kinematic  Relations 


Acceptable  nonlinear  strain-displacement  relations 
for  shells  may  be  obtained  from  one  of  two  basically  dif- 
fering approaches.   One  approach  is  to  utilize  the  finite 
strain  tensor  of  the  three-dimensional  theory  of  elasticity 
in  conjunction  with  the  Gauss-Codazzi  relations  for  a  sur- 
face [30,  31,  32,  33,  34].   Another  method  is  the  two- 
dimensional  approach  to  shell  theory  that  evades  the  ques- 
tions of  the  approximations  involved  in  the  descent  from 
three  dimensions.   This  approach  defines  strain  as  one- 
half  of  the  difference  between  the  material  and  spatial 
metric  surface  tensors.   Sanders  [35]  presented  the  first 
consistent  nonlinear  kinematic  relations  derived  by  this 
method  using  a  continuum  called  a  Cosserat  surface. 
Naghdi  [36]  discusses  both  of  the  above  methods  in  detail. 

Nonlinear  kinematic  relations  are  seldom  utilized 
in  shell  analysis  in  their  complete  form  due  to  their 
complexity.   Therefore,  simplifications  are  routinely  made. 
Throughout  the  following  analysis  it  is  assumed  that  the 
strains  in  the  middle  surface  and  that  the  rotations  about 
coordinate  axes  are  small.   These  assumptions  imply  that 
an  element  of  area  on  the  deformed  middle  surface  is  iden- 
tical in  size  to  an  element  of  area  on  the  undeformed 
middle  surface  and  that  the  difference  between  the 
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Christoffel  symbols  of  the  deformed  and  undeformed  coordi- 
nate systems  is  zero.   These  assumptions  result  in  a  con- 
siderable simplification  in  the  resulting  equations  with- 
out unduly  restricting  the  applicability  of  the  solutions. 
The  strain-displacement  equations  shown  below  were 
derived  using  the  finite  strain  description  of  the  three- 
dimensional  theory  of  elasticity  in  a  manner  similar  to 
Ogibalov  [32]  and  the  nonlinear  terms  were  regrouped  into 
surface  rotation  expressions.   The  resulting  kinematic 
rel ati  ons  are 


e^  =  (l/(l+z/R^))[U,^/A+A,j^V/AB  +  W/R^+(l/(2(l+z/R^))) 


x(-W,  /A+U/R  )^] 
^    a      a' 


eg  -  (l/(l+z/Rp))[V,p/B+B,^U/AB+W/Rg+(l/(2(l+z/Rp))) 


x(-W,^/B+V/Rj^] 


(2.3) 


^a3^  (l/(2(Uz/Rj))(V,^/A-A,3U/AB) 


+  (l/(2(l+z/Rj))(U,„/B-B,  V/AB) 


+  (l/(2(l+z/R^)(Uz/Rp)))(-W,^/A  +  U/R^)(-W,g/B  +  V/Rg). 


In  equations  (2.3)  the  covariant  components  of  the  normal 
and  shear  strains  are  denoted  by  e  ,  e„,  e  „  respectively 
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The  components  U,  V,  and  W  represent  displacements  of 

an  arbitrary  point  in  the  shell.   R   and  R„  are  the 

a       3 

principal  radii  of  curvature  and  A  and  B  are  the  coeffi- 
cients of  the  first  fundamental  form.   Comma  denotes 
partial  differentiation. 


2 . 3   Constitutive  Relations 

The  relations  between  the  components  of  strain 

and  stress  in  an  orthotropic,  linearly  elastic  material 

are  given  by  Hooke's  law  [37]  as 


e^  =  a"/E  -V  .a^/E  -V   a"/E  +a  T 
a       a   a6     3   an     n   a 


e„  =  -v^  a^'/E  +a^/E  -v^  a/E  +a„T 
6     Ba    a     B   Bn   '  n   6 


e„  =  -v„  a"/E  -V  „a^/E  +a"/E  +a  T 
n     na    a   nB    6     n   n 


2e  ,  =  a"^/G  , 

aB         at 


(2.4) 


0^    -   Bn,p 
2^Bn  -  ^   /Sn 


2e    =  a""/G 
an        an 


Young's  modulus  in  the  a  and  B  directions  is  de- 
noted by  E  and  E      respectively.   The  Poisson   Ratio  v 

a.  fi  a£ 

designates  the  contraction  in  the  a  direction  caused  by 


a  positive  normal  stress  a^  in  the  B  direction.   G  „   is 
the  shear  modulus  in  a  plane  which  is  tangent  to  the 

(a, 6)  coordinate  surface  and  (a  ,  a„)  are  the  coeffi- 

a    3 

cients  of  thermal  expansion  in  the  a  and  3  directions 
respectively.   These  relationships  are  valid  for  materials 
subjected  to  stresses  below  the  proportional  limit.   By 
introducing  Kirchhoff's  hypothesis  and  allowing  for  the 
symmetry  of  the  elastic  parameter  matrix  [38]  the  above 
equations  may  be  written  as 


^   =  t^/(^-"a3V)^^^WS^3-(^V^3V^3)T^ 

^^  =  tl/(l-v  .V   )][E  P^+v  E  e  -(a,E,  +  a  >^^  F JT]    (2.5) 


^  =  2G  ,e  „ 

aB  a3  • 


2.4   Kirchhoff's  Hypothesis 


The  reduction  of  the  three-dimensional  problem  to 
a  two-dimensional  one  requires  an  assumption  concerning 
the  variation  of  strain,  displacement,  or  stress  across 
the  shell  thickness.   To  satisfy  this  requirement  the 
fourth  assumption  in  Love's  first  approximation  [39], 
known  as  the  Kirchhoff  hypothesis,  is  introduced.   This 
hypothesis  entails  the  vanishing  of  transverse  shearing  and 
normal  strains  [33]  and  may  be  formulated  as  follows: 
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U(a,3,z)  =  u(a,3)+ze^(a,3) 
V(a,6,z)  =  v(a,6)+zep(a,3) 
W(a,3,z)  =  w(a,3) 
0  =  u/R  -w,  /A   ,   9_  =  v/Re-w,„/l 


(2.6) 


where  u,  v,  and  w  are  the  components  of  displacement  at 
the  middle  surface  in  the  a, 3,  and  normal  directions, 

respectively,  and  6   and  e„  are  the  rotations  of  the 

^       -^  a  3 

normal  to  the  middle  surface  during  deformation  about 
the  3  &  a  axes,  respectively.   The  acceptance  of  this 
assumption  is  due  to  its  great  clarity  [40].   Although 
Kirchhoff's  hypothesis  is  a  first  approximation,  its 
applicability  to  nonlinear  shell  theory  is  well  known 
[33,  35].   The  problem  of  determining  the  error  intro- 
duced by  the  hypothesis  on  the  preservation  of  the  normal 
has  thus  far  not  been  solved  exactly  [41].   Novozhilov 
[42],  Mushtari  [43],  and  Koiter  [44]  have  estimated  that 
the  errors  introduced  by  Kirchhoff's  hypothesis  are  of 
the  order  of  h/R,  which  is  small  for  thin  shells. 

Kirchhoff's  hypothesis,  modified  to  account  for 
transverse  normal  strain,  has  been  employed  in  [31,  33, 
34,  36,  45].   Effects  of  transverse  normal  strain  were 
not  included  in  the  present  analysis. 
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2  .  5   Variational  Principle 

2.5.1   Hamilton's  Principle 

The  Irish  mathematician  and  physicist,  Sir  William 
Hamilton  (1805-1865),  formulated  his  celebrated  principle 
in  dynamics  in  which  the  governing  equation  depends  ex- 
plicitly on  the  energy  of  the  system  [46].   Hamilton's 
principle  is  stated  as  an  integral  equation  in  which  the 
energy  is  integrated  over  an  interval  in  time.   In  the 
language  of  the  calculus  of  variations,  Hamilton's  prin- 
ciple states  that  the  first  variation  of  the  time  integral 
of  the  difference  between  the  kinetic  energy  T  and  the 
potential  energy  P  of  a  dynamical  system  is  zero,  that  is. 


2  'x,   % 
/  (T  -  P)dt  =  0 


(2.7) 


The  equation  is  assumed  to  hold  for  all  dynamical  systems 
whether  they  are  conservative  or  nonconservati ve .   The 
quantity  (T-P)  is  called  the  Lagrangian  function  and  is 
denoted  by  L.   With  this  designation  equation  (2.7)  be- 
comes 


/  Ldt  =  0 


(2.8) 
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and  Hamilton's  principle  then  asserts  that  the  first  vari- 
ation of  the  time-integral  of  the  Lagrangian  function  is 
zero.   Hamilton's  principle  is  employed  in  the  next  sec- 
tion to  derive  the  equations  of  motion  and  the  admissible 
boundary  conditions  for  shells. 

2.5.2   Equations  of  Motion 

Of  the  many  procedures  available  for  the  derivation 
of  the  equations  of  motion  for  a  differential  shell  element, 
Hamilton's  principle  is  superior  in  nonlinear  shell  theory 
due  to  its  efficient  treatment  of  problems  involving  curv- 
ilinear coordinates  and  because  it  gives  the  admissible 
boundary  conditions,  natural  or  forced,  that  are  to  be  used 
with  the  theory  [47,  48].   With  the  total  potential  energy 
P  expressed  as  the  difference  between  the  internal  strain 
energy  and  the  potential  energy  of  the  external  forces,  we 
obtain  from  equation  (2.8) 


5/  {/y  [(P_/2)V\  -(l/2)s''-^Y,-Jdv+/   s^  dsldt  =  0 

t,   0  ^J     s  "    ' 

I  0 


(2.9) 


Here  the  density  of  the  undeformed  body  and  the  Lagrangian 
components  of  the  displacement  vector  are  p  and  V.,  re- 
spectively.  The  symmetric  Cauchy  stress  tensor  s^"^  is 
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measured  per  unit  area  of  the  undeformed  body  and  must 
be  referred  to  the  base  vectors  in  the  deformed  body  and 
the  strain  tensor,  Y-.,  is  defined  by  equation  (2.3). 
The  components  of  the  surface  force   s^  are  referred  to 
the  base  vectors  in  the  undeformed  body  and  the  volume 
and  surface  integrals  in  equation  (2.9)  must  be  extended 
over  the  volume  and  surface  of  the  body  in  its  undeformed 
state . 

Introducing  the  kinematic,  constitutive,  and 
Kirchhoff  relations  from  previous  sections,  we  eliminate 
the  time  derivatives  of  the  variations  by  integrating  by 
parts  with  respect  to  time  and  require  the  virtual  dis- 
placements to  vanish  at  the  end  points  of  the  arbitrary 
interval  ti<t^t2,  thus  obtaining  the  equations  of  motion, 
stress  resultants  and  couples,  and  boundary  conditions. 

The  equations  of  motion  derived  in  this  section 
reflect  the  assumptions  of  Kirchhoff's  hypothesis.   Small 
strains  and  moderately  small  rotations  were  also  assumed. 
Thinness  assumptions  were  delayed  until  after  the  equa- 
tions were  synthesized. 
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The  equations  of  motion  are 


+ABp"  =  ABphu(l+h^/12R^Rg)+ABph^(l/R^+VR   )e^/12 


(BN°'^),   +(AN'^),„+B,   N^^'-A.^N^'+ABQ^/R^-ABlN^e^+N^'^e    )/R, 

LX  pLX  p  p  pCX[! 

+ABp^  =  ABphv(1+h^/12R  R   )+ABph^(1/R  +1/Rje'  /12 

UC     p  LX  p  p 


(bq"),  +(aq^),  -ab(n"/r  +N^/Rj-(Be  N^'+BeJ"^), 

ut  p  UC  p  ex,  p  ct 

-(Ae^N°'^+AegN^),g+ABp"  =  ABphw(l+h^/12R^Rg) 


(BM^) ,^+(AM^")     +A     m"^-B,^M^-ABQ" 
-  ABph^[(l/R^+1/Rp)Li+eJ/12 


(BM^^) ,^+(AM^) ,g+B,^M^"-A,gM"-ABQ^ 
=  ABph^[(l/R^+VR3)v+6p]/12 


where  ( " )  denotes  second  partial  differentiation  with 
respect  to  time  and  p"' ,  p^,  p"  represent  the  applied 
surface    1  oads . 


(2.10) 
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The  stress  resultants  and  stress  couples  are  given  by 


..a      h/2   a 

N^    -h/2  a'^'^       ^ 


,,B      h/2   3 
{%  >  =  /    {°.  }(l+z/R  )dz 
^3a     _^/2  a^" 


M°'       h/2   a 

m"^    -h/2  a"^       ^ 


M^      h/2   3 

{^«  }  -  /    {''.  }(l+z/R  )zdz 
..3a      ,  ,o   3a  ^    'a' 
M       -h/2  a 


(2.11) 


;a3 


^a3 


h/2   an   ,   ,p 

-h/2  a^"   1+z/R 

a 

h/2 

/    a«^dz 
-h/2 

h/2 

/   a^^zdz 
-h/2 


where  h  is  the  shell  thickness. 

Figure  2  illustrates  the  relations  between  the 
stress  resultants  and  couples  to  the  generalized  coordi 
nates  and  to  the  shell  element. 
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,z,w 


Figure  2.   Shell  Element 
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2.5.3   Boundary  Conditions 

From  equations  (2.9),  (2.10),  and  (2.11)  we  get 
the  following  natural  boundary  conditions  that  arise 
from  the  requirement  of  force  balance.   Along  the  edge 
of  constant  a  the  boundary  conditions  are 


N   =  N   or  u  =  u 


.,a3  ,  M^P  /  n     mO^B 

N   +M   /Rq  =  N    or  V  =  V 


Q  +M,g-N  e^-N   Bg  =  Q    or  w  =  w 


(2.12: 


M   =  M    or 


The  boundary  conditions  along  constant  3  may  be  obtained 
from  equation  (2.12)  by  interchanging  a  with  3  and  u  with 
V  . 


2  .  6   Synthesis  of  Equations 


The  usual  procedure  followed  is  to  reduce  the 
number  of  equations  and  unknowns  to  a  more  manageable 

number  by  eliminating  Q   and  Q„  from  the  equations  of 

^   a       3 

motion  (2.10)  which  are  thus  reduced  from  five  to  three 
The  force  and  moment  resultant  expressions  (2.11)  to- 
gether with  the  constitutive  relations  (2.5)  are  then 
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substituted  into  the  equations  of  motion.   Finally,  the 
strain-displacement  equations  (2.3)  and  (2,6)  are  substi- 
tuted, to  yield  three  differential  equations  of  motion 
having  u,  v,  and  w  as  dependent  variables  and  a,  3,  and 
t  as  independent  variables. 

2.6.1   Shells  of  Revolution 

In  the  preceding  analysis  R   and  R„  are  arbitrary 
but  in  all  subsequent  applications  these  principal  radii 
will  refer  only  to  shells  of  revolution.   It  can  be  shown 
[47]  that  for  a  shell  of  revolution  the  principal  radii 
may  be  expressed  as 


R„  -    -[l+(9Ro/3x3)^]^/^/(9^R^/9x3^) 


(2.13) 


h   -   %^''^'K/''3^' 


where  R  is  the  radius  of  the  latitude  circle  and  x.  is 
0  3 

the  coordinate  coincident  with  the  axis  of  revolution. 


2.6.2   Conical  Shells 


The  above  described  synthesis  may  be  carried  out 
using  equation  (2.13);  however,  the  equations  would  be 
extremely  unwieldy.   Thus,  the  procedure  will  be  followed 
only  for  a  specific  shell  of  revolution  shape  and  the  one 


28 


used  in  the  remainder  of  this  study  is  a  circular  conical 
shell.   A  conical  shell  reduces  to  a  cylindrical  shell 
when  the  semivertex  angle  becomes  zero.   Similarly,  it 
reduces  to  a  flat  circular  plate  when  the  semivertex 
angle  approaches  a  right  angle.   Therefore,  the  following 
analysis  can  handle  cones  as  well  as  the  two  related 
problems.   The  conical  coordinate  system  adopted  is  shown 
in  Figure  3.   The  coordinates,  coefficients  of  the  first 
fundamental  form,  and  principal  radii  become 


(2.14) 


In  all  subsequent  relations,  a  will  denote  the  semivertex 

angl e . 

Assumptions  relating  to  displacement  magnitudes  are  now 
introduced  into  the  analysis.   The  squares  of  inplane 
displacements  are  assumed  to  be  small  with  respect  to  the 
squares  of  normal  displacements  and  the  shell  thickness 
is  assumed  to  be  small  compared  to  the  radius. 


a    = 

s 

3   = 

H 

A    - 

1 

B    = 

s 

si  n 

a 

R      = 

a 

CO 

h- 

s 

tan 

a 
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A  -  A 


Figure  3.   Conical  Shell  Coordinates 
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Neglecting    rotary    inertia    terms,    the    nonlinear 
equations    of   motion    are 


N^+sN^  -N^+N^^  =  sphii 

s  i) 


-ctn^a  vN^/s+ctna  w,  N^/s+N?  +2N^^+sN^^+ctna  w.^N^^ 


+ctna  M?  /s+2ctna  M^^/s+ctna  M^^  =  sphv 


w,^N^+sw,g^N^+sw,gN^^-ctna  N^-ctna  v,  N^/s+w,  N^/s         (2.15) 


-etna  vN?  /s+w,   N?  /s-ctna  v,2N^^+2w,     N^^-ctna  vN^^ 


+^'/'^^'s^'^2M's+^^^'ss-f^'s^^'#/^^2M^J/s+2M^^^+sp 


=  sphw 

where 

ijj    =    e    s  i  n  a  . 
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The   nonlinear    stress    resultants    and    stress    couples 
are 


[Eh/s{l-v'^)]{su,g+s(w,^)V2+v[v,^+u+ctna  w+(w,^)V2s 


-etna  vw, ,/s]} 


^^     =  N       =   (Gh/s)(sv.  +u, ,-v+w.  w, , -etna  vw.    ) 

^  '  S         l|j  S      IJJ  s 


N       =   [Eh/s   (1-v   )]{sv,   +su+sctna  w+(w,    )   /2-ctna  vw, 

+ctn^a  v^/2+v[s^u,5+s^(w,g)^/2]}  ^^•^^' 


M^     -   [Eh^/12s(l-v^)]{-sw,^g+ctna  u ,g+ctna(w, 2)^/2 

2       2 

+v[-w,     /s-w,g+ctna(w,    )   /2s   ] 


M^^  -  M^s  =   (Gh-^/12s)(-2w,     +2w,,/s+ctna  w,  w,,/s) 
^  '^  sip         Y  s     Y 


6  3         2         2  2  2 

M       =  [Eh  /12s   (1-v   )](-w,, ,-sw,   -etna  u-ctn  a  w-vs  w,      ) 

M  S  SS 
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Substituting  equations  (2.16)  into  equations 
(2.15)  yield  three  equations  of  motion.   Noting  that 

Gh  =  K(l-v)/2 

we  have 


K{[-u/s+u,^+su,3^-Hu,^^/2s-3v,^/2s+v,^^/2-ctnaw/s 
+(w,^/2s)(w,^^-ctna  v,^)+(w,gy2s)  (w,^-ctna  v)  +  (w,^)^/2 


2   2  2    2   2   2 

-(w,^)  /2s  +ctna  vw,  /s  -ctn  a  v  /2s  +sw,gW,2^] 

+v[-u,^^/2s+v,^/2s+v,3^/2+ctnaw,^-w,^/2(w,3+w   /s 


■etna  V,  /s)+w,  /s(w,g  /2-w,  /2s-ctna  v.^+ctna  v/s) 


■etna  vw,  /2s]}  =  phsii 


K{3uy2s+u,^^/2-v/2s+v,g/2+sv,^^/2+v,^^/s+etna  w^s 
+w,^(w,^3/2+w,^^/s2)+w,^(w,^^+wys)/2 
+ctna[-sv(w,gg/2+w,  /s  )+w,g(-v-^u,,/2)+w,  u/s]/s 
+v[-u,^/2s+u,g^/2+v/2s-v,^/2-sv,^^/2+w,^w,^^/2-w,^w,^/2s 
-w,ggW,  /2+etna  u,g/s(-ctna  v+w,  )+ctna  w,^(v-u,  /2)/s 
+etna  vw,gg/2]}  =  phsv 


(2.17) 


(2.18) 
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2  2 

K{-ctna  (u+v,^+ctna  w)/s+ctna[(w,^)   /2+ww,^^^]/s  +su,gW,g^ 

+w.    (su.     +v.    ,/2+u,  ,  ,/2s)+w,  ,  ,v,,/s^+w,  ,(v,, ,/s^+v,     /2 

S  SS        SljJ  w  #     'I'  4*       #  ss 

+u,    ,/2s)+w.    ,(u,  ,/s+v.    )+u,  w,  +UW, , ,/s^+u. ,w,,/2s^ 
-vw,g^/s-v,gW,^/2s+vw,^^/2s  -v,^j^w,^/2s+ctna(-v,^  /s   -uv,   /s 

-vV'^/s^-U'//2s^-v,^2/2-u,^v,^/2sH-vv,^/s-vv,^^/2-u,^^v/2s 

2       2  2  2 

-V  /2s   )+ctn  a  (-v,,w/s   )+v[-ctna  u,  +w,      (v,  +ctna  w) 

^^>(^'s^/2s-v,3s/2)+u,3W,^^/s+w,^(v,^^/2+ctnaw,^/2-u,^^/2s) 

+uw,^^+u,^w,^-v,^w,^^-u,^w,^^/s+vw,^^/s+u,^w,^/2s2  (2.19) 

+v,^w,^/2s-vw,^/2s2+v,^w,^/2s+ctna  (-u,^v,^/s-u,^^v/2s 

+v,^/2+u,^v,2/2s+vv,2^/2-u,^v/2s^-vv,^/s+v^/2s^)]} 

"°^-^^'ssss-2^'ss#/^-^w/^'^2w,^^/s2-2w,333-4w,^^/s3 

+v[(ctna/s^)((-7/2s)(w,^w,^^)+3(w,^)^/s^+w,2W,^^/2s 
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where  the  membrane  stiffness  is 


K  =  Eh/d-v'^) 


(2.20) 


and  the  flexural  stiffness  is 


D  =  Eh^/[12(l-v^)] . 


(2.21) 


CHAPTER  3 

NONLINEAR  DYNAMIC  STABILITY 
EQUATIONS  OF  CONICAL  SHELLS 

The  nonlinear  dynamic  stability  equations  are 
obtained  from  the  nonlinear  equations  of  motion  by 
employing  a  perturbation  analysis.   Substituting  the 
displacements 


u  =  Uy  +  u. 


V  =  Vj  +  V( 


(3.1) 


w  =  w.  +  w. 


and  the  stress  resultants  and  stress  couple  expressions 


1^  =  N^  +  Nj 


M^  =  Mj  +  Mg 


M''  =  Mj  +  Mj 


(3.2) 


N^^  =  Nf  +  Nf 


M^Q  ^  Mf  +  Mf 

i       D 


35 


36 


and  the  load  relation 


P  =  Pj  +  P[ 


(3.3) 


into  equations  (2.15)  results  in  two  coupled  sets  of  equa- 
tions.  The  quantities  with  subscript  I  are  measured  from 
the  undeformed  state  and  the  quantities  with  subscript  B 
are  small  but  finite  perturbations.   This  gives  one  equa- 
tion set  describing  the  motion  in  the  stable  space  surround 
ing  the  undeformed  state  and  another  set,  the  stability 
equations,  describing  the  motion  during  passage  thru  the 
nearest  bifurcation  point. 

The  resulting  set  of  nonlinear  dynamic  stability 
equations  are  as  follows  where  the  subscript  B  has  been 
dropped  for  convenience: 


K[(-u/s+u,2+su,^2+u,^/2s-3v,^/2s+v,g^/2-ctna  w/s] 


+^(-^'#/^s+v,^/2s+v,^^/2+ctna  w,^) 
^(^'s^V^^+^'s/V^^-^^'s/^-^'^/^^^-^^'^'s'^'ss^ 


(3.4) 


+v(-w,^/2-w,^w,^^/2s+w.^w,^^/2s-w,^/2s^)]  =  phsu 
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+^(-^V2s-.^'s/2+^/2s-v,3/2-sv,3^/2) 

+  (w,^w.^3/2+w,^w,^^/s^w,^w,^^/2+w,^w,^/2s+ctn^aww,^/s^) 


+v(w,^w.    ,/2-w,  w,,/2s-w,     w,,/2)]  =  phsv 

^      S      Slj^  s      tl^ SS      Y 


K{-ctna(u+v, ,+ctna  w)/s+v(-ctna  u,    ) 
+ctna  [w,^/2s^+ww,^^^/s^+v(w,5/2+ww,gg)]} 


-4w.M/s^+w,^^/s-w,^/s2)+ctna[w,J+w,^3W       /s^ 


2     2  3  3  2 

+w,^,/s  +w,^w.  ^^-w, ,w.    ,/s  -w.  w,,,/s  +w,,w.  ^,/s 


+w.  w,    ,  ,/s  +v(-3w,,w,^,/s  +3w,,/2s  +w.  w,,  ./s"^ 


■^'ss^  V  -^'s^'s^/'   )^^  ^  ^P" 


+w,  N^+sw.   J]t+sw,  Nt,  +w,, ,Nt/s+w,,Nt,,/s  =  phsw 

S    I  SS    I  S    Is        'd)'>iJ    I  \h    I    Uj 


(3.5) 


(3.6; 
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The  single  underlined  quantities  represent  non- 
linear membrane  stiffness  terms.   The  double  underlined 
quantities  represent  nonlinear  bending  stiffness  terms. 
The  nonlinear  stability  terms  are  denoted  by  a  dashed 
underline.   These  equations  embody  the  assumptions  that 


and  that 


Wy   <  Wr 


Un  <  Wr 


Vn  <  w. 


In  the  present  study,  a  thin  truncated  conical 
shell  is  loaded  by  a  constant  axial  load  and  a  nearly 
axisymmetric  Heaviside  pressure  load.   The  membrane- 
state  stress  resultants  are  given  by 

Nj  =  -P/(us  sin2a)-pj(s tana) (l-cosAQt)/2 

fi      n 
Nj  =  -Pj(stana)(l-cosA  t) 


where  P  is  the  applied  axial  load,  p"  is  the  magnitude 
of  the  Heaviside  pressure  load,  and  A   is  the  natural 
frequency  of  the  breathing  mode  of  the  shell. 


CHAPTER  4 
METHOD  OF  SOLUTION 

4  . 1   Circumferential  Modal  Analysis 

The  nonlinear  partial  differential  equations 
derived  in  previous  chapters  are  spatially  two-dimen- 
sional.  The  two  directions  are  meridional  and  circum- 
ferential.  In  this  analysis,  the  two-dimensional  equa- 
tions are  reduced  to  one-dimensional  equations  by 
utilizing  circumferential  modal  analysis.   The  inde- 
pendent variable  6  is  eliminated  by  expanding  all 
dependent  variables  into  sine  or  cosine  series  in  the 
circumferential  direction.   As  a  result  of  the  trig- 
nometic  series  expansions,  there  is  one  set  of  govern- 
ing equations  for  each  circumferential  wave  number 
considered.   For  a  linearized  analysis  the  sets  are  un- 
coupled.  For  a  nonlinear  analysis  the  equation  sets  are 
modally  coupled  through  the  quadratic  terms.   This  cou- 
pling arises  when  use  is  made  of  the  trignometric  prod- 
uct identities 


cos(ie)cos( j0)  =  (l/2)[cos(i-j )e  +  cos(i+j  )0] 
sin(i6)sin(je)  =  ( 1 /2 ) [cos ( i - j )e-cos ( i+j )e] 


(4.1) 
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in  order  to  facilitate  the  equating  of  like  coefficients 
of  Fourier  expansion  terms  which  is  required  to  identify 
the    equation    sets. 

Substituting    the    expressions 


u(s,i|j,t)         ==       E        u    (s,t)    cos(n;|;/sina) 
n  =  o        "^ 
Q 
v(s,i|i,t)         =      E        v(s,t)    sin(ni|j/sina) 
n  =  o 
Q 
w(s,i|j,t)         =      T.        w    (s,t)    cos(nifj/sina) 
n  =  o         " 
Q 
p(s,4^,t)         =      E        p    (s,t)    cos  (ni]j/sina) 
n  =  o         '^ 


(4.2) 


into  the  nonlinear  dynamic  stability  equations  (3.4  ■ 
3.6)  and  equating  like  coefficients  of  Fourier  expan^ 
sion  terms  yield  the  following: 


(-l/s){l+[n^(l-v)/(2sin^a)]}u  +u  ,  +su  , 


+(n/2s  sina) (-3+v)v^+(n/2sina) (l+v)v^,^-(ctna  w^/s)  ,^^^ 


+vctna  w  ,   =  phsii  /K 
n  s    ^    n 
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(n/2s  sina)(-3+v)u  - ( n/2si na) ( 1 +v ) u  , 


+{-[(l-v)/2]-(n2/sin2a)}(v^/s)+[(1-v)v^,^/2] 


(4.4) 


+[(l-v)sv  ,   /2]-[n  etna  w  /(s  sina)]  =  phsv  /K 


-etna  [(u^/s)+vu^  ,^  +  (n/sina) (v^/s)  +  {ctna  w^/s)] 

+  (D/K){[(4n^/sin^a)-(n^/sin^a)](w^/s-^)  +  [(2n^/sin^a)+l] 

x[(w    ,      /s)-(w    ,    /s^)]-2w    ,         -sw    ,  }+sp    /K 

^    n    ss  n    s         ' -'         n    sss        n    ssss        ^n 


-[P.w    ,      /(K7Tsin2a)]  +  (1 -cosX    t)[-s    tanc 
"-Inss  '  -I    \  0 


c    s 


(s^tana6",g5/2)+(2n^6"/sin2a)-(2n3"/sina)]/K 


+n^{ctna{(-n^w^/4s^sin^a)+v[(w^,^/2)+(w^,^/4) 


(4.5) 


^Vo'ss^(Vn'ss/2):i>-^(h2ctna/12){w^,jH-(w^,J/2) 

0  0         0  o  o 

-[n    /(2s  sin  a) ] (w^w^ ,^^-w^ , ^ ) } }+n^{ctna[ ( -n    w^w^/ 

/    2    .    2    ^  X 
(  s    sin    a) ) 


+v(w„,  w  ,  +w  ,   w  +w  w  ,   )]+(h  ctna/12)[2w  ,   w  , 
OS  n  s   0  ss  n   o  n  ss  -"  ^      '   "-   o  ss  n  ss 


2   2    2 
(n  /s  sin  a)w  ,   w  ]}  =  phsw  /K 
^  '  0  ss  n-^        n 
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where 


n    =    circumferential    wave    number 


1    n  =  o 
0    n>o 


I    w./2[n    p  /  .  ,     x+i_i    n,.        i] 
r     "-      ^(  1+n)         ^    1  -n  M 


,   n 


Z    w./2[n    D/.^     x+y^|i-n|p|.       ,] 


0  n  =  o 

1  n>o 


(4.6) 


,1    n>o 


1    i>n+l 
y     =  "^  0    i  =  n 
■1    im-l 


4.2   Elimination  of  Meridional  Derivatives 


The  next  phase  of  the  analysis  involves  the  re- 
duction of  the  partial  differential  equations  sets  result- 
ing from  the  circumferential  modal  analysis  to  ordinary 
temporal  differential  equation  sets  which  are  then  amen- 
able to  numerical  solution. 

Referring  again  to  Figure  1  it  is  seen  that  the 
available  solution  methods  for  this  task  fall  generally 
into  four  categories.   These  are  finite  element  methods, 
finite  difference  methods,  variational  principles,  and 
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methods  of  weighted  residuals.   As  noted  by  Finlayson 
[49],  for  certain  types  of  linear  problems  these  methods 
can  be  shown  to  be  equivalent  to  each  other.   From  the 
literature  abounding  on  shell  dynamics  one  can  safely 
conclude  that  any  of  these  methods  when  judiciously  ap- 
plied will  give  acceptable  results.   There  are,  however, 
marked  differences  between  the  methods  in  the  amount  and 
quality  of  manual  mathematics  and  the  quantity  of  machine 
computing  time  required  to  give  equivalent  results. 

The  desired  characteristics  of  the  solution  method 
for  the  particular  problem  of  interest  in  this  disserta- 
tion are  now  discussed.   The  solution  method  utilized 
should  yield  an  approximate  direct  solution  where  con- 
vergence would  be  monotonic  and  assured.   Spatial  control 
of  error  residuals  was  desired.   Because  the  taper  in  the 
conical  shell  results  in  n on  symmetric  stiffness  matrices, 
a  solution  technique  applicable  to  non-self-adjoint  prob- 
lems was  required.   Finally,  it  was  hoped  that  the  chosen 
solution  method  would  result  in  short  computer  compile 
times  because  experience  has  shown  that  the  majority  of 
computer  time  required  for  such  a  problem  is  in  successive 
compile  times  during  development  and  debugging. 

In  matching  these  desired  characteristics  to  the 
known  characteristics  of  the  classical  solution  methods 
an  elimination  process  was  begun  which  finally  resulted 
in  the  algorithm  used  herein. 
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Finite  element  techniques  were  eliminated  from 
consideration  because  this  method  is  not  competitive  with 
other  methods  in  problems  with  regular  geometries  and 
formulatable  boundary  conditions  lying  along  coordinate 
lines.   This  method  requires  significantly  more  manual 
mathematics  and  computer  execution  time  than  equilibrium 
equation  solution  methods  in  order  to  give  comparable 
resul ts  . 

Finite  difference  methods  have  been  shown  to  be 
extremely  powerful  in  solving  nonlinear  dynamic  stability 
problems  involving  shells.   This  method  has  been  used  with 
particular  success  when  coupled  with  circumferential  mode 
analysis.   It  was  shown,  however,  by  Radosta  [50]  that 
convergence  is  slow  when  a  meridional  mesh  is  used  with 
a  shell  under  high  axial  load.   In  the  particular  problem 
at  hand,  i.e.  a  conical  frustum  subject  to  a  high  axial 
load  combined  with  a  time  dependent  pressure  load,  it  was 
felt  that  a  suitable  selection  of  number  and  spacing  of 
meridional  mesh  points  would  degrade  into  a  highly  time 
consuming  trial  and  error  effort. 

The  final  phase  of  solution  method  selection  con- 
sisted of  choosing  one  of  the  four  methods  of  weighted 
residuals:   Galerkin,  Least  Squares,  Collocation,  and  Sub- 
domain.   The  least  squares  method  was  not  seriously  con- 
sidered because  it  yields  ordinary  differential  equations 
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of  the  second  degree  in  the  temporal  operator.   The 
collocation  technique  was  eliminated  because  it  enhances 
the  same  point  location  selection  drawbacks  as  does  the 
finite  difference  method.   Again,  this  is  normally  an 
insignificant  problem  but  was  made  significant  here  by 
the  expected  buckled  states  resulting  from  the  combined 
type  loadings  being  considered.   The  final  choice  between 
the  Galerkin  and  subdomain  methods  was  based  on  the  method 
of  error  residual  control.   In  the  Galerkin  technique  the 
error  produced  by  each  individual  term  in  the  assumed  solu- 
tion is  forced  to  average  to  zero  over  the  entire  length 
of  the  shell.   In  the  subdomain  method  the  error  in  the 
total  assumed  solution  is  forced  to  average  to  zero  over 
each  subdomain.   This  type  of  spatial  control  for  error 
residual  was  considered  superior  to  the  modal  type  of 
error  control  offered  by  the  Galerkin  method  and  there- 
fore, the  subdomain  technique  was  selected  for  use  in  this 
investigation. 

4.2.1   Subdomain  Method 

A  description  of  the  theoretical  basis  of  the  sub- 
domain  method  is  now  in  order.   Consider  the  equation 


L[u(x,t)]-0    ,inV  (4.7) 

where  L  is  a  nonlinear,  non-sel f -adjoi nt ,  ordinary 
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differential  operator.   If  the  approximation  solution  is 
represented  in  the  form 


Z  a.(t)(l).(x) 
i  =  l  ^     ^ 


(4.8) 


and  this  is  substituted  into  equation  (4.7),  we  obtain 
an  approximation 


L[u(x,t)]  =  e 


(4.9) 


where  e  is  the  error.   We  now  impose  constraints  on  e  by 
requiring  it  to  average  to  zero  over  a  subdomain  v  of 
the  total  domain  V 


/L[u(x,t)]dv  =  0  . 
v 


(4.10) 


The  number  of  subdomains  v  are  chosen  to  equal  m  the 
number  of  time  dependent  unknowns  in  (4.8).   The  differen- 
tial equation,  integrated  over  the  subdomain   is  then  zero 
hence  the  name  subdomain  method.   As  m  increases,  the  dif- 
ferential equation  is  satisfied  on  the  average  in  smaller 
and  smaller  subdomains,  and  presumably  approaches  zero 
everywhere.   This  process  yields  a  dynamically  coupled 
second  order  set  of  nonlinear  temporal  differential  equa- 
tions which  may  in  turn  be  solved  by  numerical  methods. 
Equation  (4.8)  is  called  the  trial  function  and 
may  be  selected  in  a  variety  of  ways.   It  may  be  chosen 
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to  satisfy  the  boundary  conditions,  but  not  the  differen- 
tial equations  or  to  satisfy  the  differential  equations 
but  not  the  boundary  conditions.   In  order  to  simplify 
the  required  hand  calculations,  a  set  of  complete  power 
series  were  chosen  as  the  trial  functions  for  this  prob- 
lem.  The  trial  functions  were  then  forced  to  satisfy 
the  boundary  conditions  by  a  modified  application  of  the 
Lagrangian  multiplier  method  which  can  be  more  properly 
thought  of  as  a  coordinate  transformation  to  require  an 
assumed  power  series  solution  to  satisfy  the  essential 
boundary  conditions.   An  approximate  satisfaction  of  the 
differential  equations  was  then  obtained  by  the  subdomain 
method . 

4.2.2  Method  of  Satisfying  Boundary  Conditions 

A  description  of  the  technique  utilized  to  satisfy 
boundary  conditions  will  now  be  presented.   Consider  the 
trial  function  (4.8)  in  matrix  form. 


u  =  L  cj)J  {a} 
(lxm)(mxl) 


(4.11) 


and  the  j  essential  boundary  constraints, 


[Cj]{a}  =  0 
(jxm)(mxl) 


(4.12) 
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Equation  (4.12)  may  be  partitioned  as  follows 


[     c^      I    4  ]  r   q^ 


[jx(m-j)]   (jxj)   1  q 


■(mxl) 


0. 


(4.13) 


Equation  (4.13)  may  be  rearranged  as 


{q2>  =  -iC^rh       C^    ]   {q^} 
(jxl)     (jxj)  [jx(m-j)][(m-j)xl] 


(4.14) 


and  noting  that 


{q^}    =    [iKq^} 


(4.15) 


we    obtain 


{a}      =    [        B      ]      {q^} 
(mxl)         [mx(m-j)][(m-j)xl  ] 


(4.16) 


and    (4.11)    becomes 


u    =    L(J)J[B]{q^} 


(4.17) 


where  [B]  may  be  considered  as  a  projection  operator 
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After  nondimensionalizing  equations  (4.3  -  4.5) 
according  to 


2  n 

s„v 
2  n 

s„w 
2  n 


T  =  /E/p(l-v2)  t/s, 


p^  =  Eh/(l-v^)s2  P^ 


X^    =  S2/p(l-v2)/E  A^ 


=  S2h 
I     cr  I 


P    =  2TTEh  cos  a//3(l-v2 
cr  ^ 


(4.18) 


derivatives  are  likewise 
nond i mens ionali zed. 


the  trial  functions  are  introduced  into  the  equations 
The  trial  functions  are 


«•                  - 

~\ 

\ 

\ 

Pn 

V.                 J 

N-1 
}    =      ^      ( 
i  =0 


c  . 
m 


m 


>    s 


(4.19) 


After  applying  the  transforms  satisfying  the 
boundary  conditions,  the  equations  are  then  integrated 
over  (N-j)  subdomains  and  the  average  residual  error  for 
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each  subdomain  equated  to  zero.   During  the  original 
formulation  it  was  noted  that  this  approach  leads  to  a 
singular  mass  matrix  for  equally  spaced  subdomains.   To 
circumvent  this  problem  the  equations  were  mutliplied 
through  by  s.   This  is  similar  to  a  procedure  known  as 
the  method  of  moments  first  employed  by  Polhausen  [51]. 

Having  completed  the  above  sequence  of  operations 
results  in  Q  sets  of  nonlinear  second  order  temporal  dif- 
ferential equations.   Each  set  contains  3  (N-j)  equations 
and  the  sets  are  coupled  in  the  nonlinear  terms.   In  ma- 
trix form  a  typical  set  may  be  expressed  as 

[m][B]{q}^"^!:k][B]{q}("^[D][B]{q}^"^{NI.}^"^  -  {P}^"^ 

(4.20) 


where 


[m]  -  mass  matrix 

[B]  -  matrix  requiring  power  series  solution  to 
satisfy  the  prescribed  boundary  conditions 


{q} 


n) 


generalized  coordinates  of  the  nth  cir- 
cumferential mode 

[k]  -  stiffness  matrix 

[D]  -  stability  matrix 

NL.  -  LqJ[B  ]  [C.][B  ]{q  }-  power  series  products 


{P} 


load  vector 
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As  an  example,  for  a  problem  retaining  two  circumferen- 
tial modes  with  eight  subdomains  along  the  meridian, 
equations  (4.20)  would  represent  a  total  set  of  forty- 
eight  coupled  nonlinear  second  order  ordinary  differen- 
tial equations  to  solve. 


CHAPTER  5 
NUMERICAL  ANALYSIS  FOR  CONICAL  SHELLS 

The  sets  of  equations  (4.20)  were  coded  into 
computer  language  for  numerical  solution.   A  particular 
shell  specimen  was  selected  for  study.   The  salient 
geometry  and  boundary  conditions  for  the  shell  studied 
are  given  in  Appendix  I.   Fortran  IV,  level  H  machine 
language  was  used  and  the  program  was  constructed  for 
double  precision  arithmetic.   A  flow  diagram  of  the  oper- 
ational characteristics  is  given  in  Appendix  III.   A  list- 
ing of  the  program  is  given  in  Appendix  IV. 

5. 1   Functional  Description  of  the  Computer  Program 

If  inplane  inertia  is  considered  ignorable  with 
respect  to  normal  inertia,  then  equations  (4.20)  separate 
into  coupled  sets  of  2(N-j)  nonlinear  algebraic  equations 
and  2(N-j)  nonlinear  first  order  differential  equations. 

The  flow  of  computations  proceeds  as  follows.   The 
linear  and  nonlinear  stiffness  matrices,  mass  matrices, 
boundary  transform  matrices,  buckling  matrices,  and  load- 
ing vectors  are  generated  by  the  program.   At  time  equal 
to  zero  the  problem  solution  begins  with  known  external 
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loading  and  initial  conditions.   The  known  normal  deflec- 
tions are  substituted  into  the  first  two  equations  of 
equilibrium  and  the  resulting  u,  v  vectors  are  generated 
utilizing  a  Gauss  elimination  technique  with  complete 
pivoting.   If  nonlinear  effects  in  the  first  two  equa- 
tions of  equilibrium  are  considered,  this  mode  is  itera- 
tive in  nature.   Having  generated  the  u's  and  v's  for 
time  zero,  they  are  substituted  into  the  sets  of  3rd 
dynamic  equilibrium  equations  for  temporal  integration. 
The  time  integration  is  accomplished  by  a  Hamming    pre- 
dictor-corrector integration  technique.   This  generates 
the  w's  for  the  next  time  step.   Normal  velocities  are 
also  computed.   This  completes  the  sequence  of  advance 
to  the  first  time  increment  beyond  time  zero  and  provides 
the  required  initial  conditions  in  order  to  begin  the 
next  advance  in  time.   The  sequence  then  recycles  and  the 
temporal  advance  continues.   At  any  time  interval  the  an- 
alyst chooses  the  solved  for  dependent  variables  are 
transformed  from  the  solution  space  to  real  space  and 
displayed  as  output.   Each  output  consists  of  the  com- 
plete meridional  profile  for  each  mode  in  the  solution. 

The  build  up  of  errors  in  the  numerical  solution 
is  controlled  by  the  following  mechanisms.   A  loss  of  sig- 
nificance control  indicator  is  incorporated  into  the  Gauss 
elimination  routine  which  warns  the  analyst  when  solutions 


54 


are  being  derived  from  a  nearly  singular  coefficient 
matrix.   An  error  bound  control  is  available  in  the 
predictor-corrector  integrator  which  limits  this  type 
of  error  by  automatically  adjusting  the  size  of  the 
time  step.   Round  off  errors  are  minimized  by  making 
all  calculations  accurate  to  sixteen  significant  figures 
Such  rigorous  control  of  possible  numerical  errors  is 
essential  in  dynamic  stability  analysis  because  of  the 
possibility  that  numerical  instability  might  be  mistaken 
for  actual  structural  instability. 

5 .  2   Checking  the  Computer  Program 


The  control  of  clerical  and  conceptual  errors  in 
the  construction  of  the  computer  program  was  done  by  con- 
tinually checking  back  to  simpler  known  cases.   In  some 
instances  the  known  cases  could  be  obtained  by  hand  cal- 
culations, in  others,  the  results  of  other  analysts  were 
used.   For  the  shell  specimen  studied,  the  computer  pro- 
gram has  been  verified  by  checking  against  known  solutions 
for  static  stability,  static  deflection,  free  vibration 
modes  and  frequencies  [39],  linear  dynamic  response,  and 
regions  of  parametric  instability  [52].   A  convergence 
study  was  also  conducted  to  determine  the  convergence 
properties  of  the  solution  for  increasing  numbers  of 
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subdomai ns .   The  results  of  this  study  are  presented  in 
Appendi  x  II. 

5  .  3   Dynamic  Stability  Results 

After  the  above  verifications  were  completed,  the 
program  was  used  to  study  the  linearized  and  the  nonlinear 
dynamic  stability  of  a  short  conical  frustum  with  simply 
supported  boundaries  subject  to  various  levels  of  constant 
axial  loads  and  critical  levels  of  nearly  axisymmetric 
Heaviside  pressure  blasts.   Eight  subdomains  along  the 
meridional  direction  were  used.   The  driven  breathing  mode 
and  an  arbitrary  preferred  flexural  mode  were  retained  in 
the  analysis.   Complete  coupling  between  the  modes  was  re- 
tained through  the  quadratic  terms  for  the  nonlinear  case. 

The  critical  static  normal  pressure  for  the  speci- 
men studied  is  42.8  pounds  per  square  inch  and  the  critical 
circumferential  wave  number  is  eight.   This  represents  a 
mid-span  deflection  in  the  breathing  mode  of   9   percent 
of  the  thickness.   The  ordinate  scale  (W)  of  all  figures 
depicting  time  history  responses  is  normalized  by  this 
critical  deflection.   The  critical  static  axial  load  is 
40,400  pounds.   The  applied  loads  for  time  history  graphs 
will  be  noted  as 

(a,  b,  c) 
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where 


a  =  P/P   ,  b  =  p  /p   ,  c  =  p  /p 


Additionally,  note  that 


where 


p   =  p  H(t' 

•^0     ^0 


p   =  p  H(t 
^n    '  n  ^ 


H(t) 


denotes  a  Heaviside  pressure  pulse 


5.3.1   Dynamic  Modal  Prefer e n c e 


The  program  was  used  to  determine  the  critical 
flexural  modes  during  dynamic  instability.   The  results 
are  shown  in  Figure  4  and  the  critical  static  buckling 
modes  are  indicated.   Because  only  a  few  closely  spaced 
critical  modes  dominate  the  flexural  response  and  because 
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coupling  between  these  modes  is  very  weak  [10],  the  param- 
eter study  that  follows  is  reduced  by  considering  the  in- 
stability of  only  a  sequence  of  single  flexural  modes  to 
find  the  mode  indicating  instability  at  the  lowest  load. 

5.3.2   Linearized  Dynamic  Stability 

The  linearized  dynamic  analysis  was  obtained  by 
considering  the  vector  {NL}^"^  in  equation  (4.20)  to  be 
zero.   This  means  that  the  internal  loads  from  the  pre- 
buckling  state  may  interact  with  postbuckling  deflections 
to  cause  primary  dynamic  instability  (i.e.  snap  through) 
but  internal  parametric  instabilities  due  to  the  pulsat- 
ing nature  of  the  membrane  state  are  prohibited. 

The  critical  circumferential  wave  number  for 
dynamic  instability  was  found  utilizing  this  analysis  and 
the  results  are  shown  in  Figure  5. 

While  conducting  a  series  of  linearized  dynamic 
stability  computer  runs  at  various  combinations  of  external 
loadings,  a  feature  of  the  numerical  model  was  discovered 
which  represents  an  improvement  over  other  computer  programs 
used  for  this  purpose.   This  feature  is  that  the  meridional 
profile  undergoes  a  marked  and  rapid  change  of  character 
during  runs  at  or  above  the  critical  dynamic  load.   For 
loadings  below  critical,  the  meridional  profile  pulsates  at 
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one  half  wave,  slightly  skewed  due  to  the  conical  taper. 
At  or  above  the  critical  load  the  meridional  profile 
begins  to  pulsate  as  a  half  wave  but  later  snaps  to  two 
or  more  half  waves.   The  time  of  snap,  number  and  shape 
of  the  waves  in  the  buckled  meridional  profile  depend  on 
the  axial  load  to  normal  pressure  ratio  and  their  rela- 
tion to  static  critical  loads.   Examples  of  this  are  shown 
in  Figures  6  and  7.   It  is  felt  that  this  improved  ability 
of  detecting  dynamic  instabilities  results  from  two  areas 
of  the  analysis.   One  being  the  spatial  error  residual  con- 
trol offered  by  the  modified  subdomain  method  employed  and 
the  other  being  the  strict  control  of  error  bound  in  the 
numerical  work.   Additionally,  it  should  be  noted  that 
most  published  works  on  this  type  of  problem  utilize  one 
term  assumed  solutions  for  preferred  modes  and  clearly 
could  not  see  this  effect  in  their  analysis. 

A  large  number  of  linearized  dynamic  stability 
runs  were  made  at  various  combinations  of  static  axial 
and  Heaviside  pressure  loads.   The  results  were  collated 
into  a  dynamic  stability  interaction  curve  and  are  shown 
in  Figure  10. 

5.3.3   Nonlinear  Dynamic  Stability 


In  the  case  of  nonlinear  dynamic  stability,  the 
possibility  of  autoparametri c  excitation  of  preferred 
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flexural  modes  by  the  pulsating  breathing  mode  is  allowed, 
in  addition  to  snap  through  or  primary  dynamic  instability. 
The  question  here  is  whether  such  excitation  will  occur 
during  the  primary  instability  motion  of  the  preferred  mode 
sufficiently  to  affect  the  magnitude  of  the  external  loads 
required  to  induce  dynamic  instability.   In  contrast  to  the 
external  parametric  excitation  problems,  the  total  energy 
of  the  shell  for  the  present  problem  is  constant  for  all 
t>0.   Thus,  the  unstable  modes  when  growing  must  derive 
their  energy  from  the  breathing  mode  with  a  corresponding 
decrease  of  energy  in  that  mode.   To  evaluate  this  inter- 
play a  full  nonlinear  analysis  is  necessary. 

Conceptually,  the  third  equation  of  nonlinear  dy- 
namic stability  may  be  expressed  in  an  average  sense  as 
fol 1 ows : 


"o  '   K^o   *   ^<%'  *   f(«n)  =  Po«(') 


w,  +  [A^  -  g(wj]w^  +  1(^'%)  +  h(w^)w^ 


(5.1  a,b) 


P,H(t) 


where  u^,  w   represent  breathing  mode  deflections  measured 

0     0 

from  the  undeformed  state  and  u  ,  v  ,  w   represent  the 

n   n   n    ^ 

poststable  growth  of  the  buckling  mode.   A   and  A   denote 
the  linear  eigenvalues  of  the  respective  modes  and  f,  g,  h. 
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k,  and  1  are  differential  operators.   The  superscript  '^ 
denotes  an  average  over  the  shell  length.   Since  w   is 

3  ^0 

excited  by  a  Heaviside  pulse,  the  w   response  is  initially 
periodic  and  equation  (5.1  b)  takes  the  character  of  an 
i nhomogeneous  Mathieu  equation.   The  primary  buckling 
terms  denoted  by  a  single  underline  cause  an  insignificant 
shift  to  the  left  on  a  parametric  stability  diagram.   In 
the  absence  of  imperfections  in  the  preferred  flexural  mode, 
denoted  by  the  double  underlined  term,  equation  (5.1  b) 
lies  well  within  a  stable  zone  of  a  Mathieu  diagram  and  no 
significant  parametric  growth  of  the  buckling  mode  would 
be  expected. 

For  finite  deterministic  imperfections,  however, 
the  location  of  stability  boundaries  becomes  more  difficult 
and  for  significantly  large  p   H(t),  solutions  to  equation 
(5.1  b) may  exhibit  parametric  growth  at  or  below  the  dynamic 
loads  required  to  produce  immediate  snap  through  and  may 
lead  to  a  delayed  dynamic  snap  through.   Whether  or  not 
this  beating  type  of  growth  occurs  with  sufficient  quick- 
ness to  have  an  effect  on  the  primary  dynamic  instability 
can  be  determined  only  by  direct  temporal  integration  of 
two  coupled  sets  of  equation  (4.20). 

A  series  of  runs  with  the  computer  program  was  made 
to  determine  the  significance  of  this  effect.   Typical  re- 
sults are  shown  in  Figures  8  and  9.   As  in  the  case  of 
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linearized  dynamic  stability,  a  large  number  of  computer 
runs  were  made  at  various  combinations  of  static  axial  and 
Heaviside  pressure  loads  and  the  results  collated  into 
a  nonlinear  dynamic  stability  interaction  curve  shown 
in  Figure  10. 

The  linearized  dynamic  stability  curve  lies  below 
the  static  stability  curve  in  Figure  10  due  to  the  in- 
creased severity  of  a  suddenly  applied  load  over  that  of 
a  statically  applied  load.   The  peak  membrane  forces  de- 
veloped by  the  breathing  mode  which  is  excited  by  a  Heavi- 
side pulse  are  twice  as  large  as  that  of  an  equivalent 
static  case  and  contribute  more  to  the  developemnt  of 
instability  in  preferred  flexural  modes  than  a  static 
situation.   This  is  known  as  the  dynamic  load  factor  effect 
The  nonlinear  dynamic  stability  curve  lies  below  the 
linearized  dynamic  stability  curve  because  internal  vibra- 
tions interact  between  coupled  modes  so  as  to  produce  un- 
stable beating  resonances  in  the  preferred  buckling  mode. 
Such  unstable  vibrations  lead  to  a  delayed  dynamic  snap 
through  and  this  phenomenon  is  called  autoparametric  ex- 
ci  tati  on . 
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Figure  10.   Dynamic  Stability  Interaction  Curves 


CHAPTER  6 
RESULTS  AND  CONCLUSIONS 

In  the  present  analysis,  a  nonlinear  shell  theory 
is  derived  and  employed  to  study  the  nature  of  nonlinear 
dynamic  instability  of  a  truncated  thin  circular  conical 
shell  structure  which  is  considered  to  be  loaded  by  a 
constant  axial  load  and  a  nearly  a xi symmetric  Heaviside 
pressure  load.   A  solution  method  is  developed  which  sat- 
isfies a  boundary  condition  exactly  and  which  converges 
toward  the  exact  solution  of  the  governing  equations  with 
increasing  subdomains.   The  method  avoids  the  necessity 
of  assuming  the  shapes  of  prebuckling  or  postbuckling 
meridional  profiles. 

The  results  of  this  study  confirm  the  feasibility 
of  the  method  of  solution  developed  for  this  analysis.   The 
method  converges  rapidly  yet  can  be  applied  to  nonlinear 
problems  with  minimum  amounts  of  manual  mathematics.   The 
trial  functions  are  simple  to  manipulate  yet  satisfy  any 
formulatable  boundary  condition.   The  method  is  particu- 
larly suitable  for  problems  where  poststable  modes  are  not 
known  in  advance  as  is  usually  the  case  for  combined  dynamic 
loadings. 
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For  the  particular  loading  conditions  studied, 
the  results  indicate  that  dynamic  instability  will  occur 
at  loads  below  critical  static  loads.   This  reduction  in 
critical  dynamic  loads  was  shown  to  be  the  result  of  both 
the  dynamic  load  factor  effect  and  of  autoparametric  ex- 
citation.  The  interaction  curves  which  are  generated 
illustrate  these  effects  quantitatively  under  conditions 
of  combined  dynamic  loading. 

A  criterion  for  dynamic  buckling  is  established 
based  on  meridional  mode  shape  changes.   This  ability  to 
detect  sudden  jumps  in  the  meridional  profile  aids  to 
verify  instability  detected  by  divergence  on  displacement 
time  history  curves  and  provides  additional  information 
about  the  poststable  state. 


APPENDIX    I 


APPENDIX  I 
SHELL  SPECIMEN  AND  BOUNDARY  CONDITIONS 

The  conical  frustum  utilized  in  this  investigation 
has  the  following  properties  (see  Figure  3): 

Material  -  1020  steel 

a  =  20° 

s-j  =  5.85  inches 

S2  =  14.22  inches 

R-l  =  2.0   inches 

R2  =  4.863  inches 

h  =  0.02  inches 

V  =  0.3  . 

The  boundary  conditions  assumed  for  the  analysis 


u  =  v  =  w  =  M   =0   ats-,  and  Sr,   . 


This  results  in  eight  constraint  equations  which  when  ex- 
pressed in  terms  of  power  series  solutions  for  eight  sub- 
domains  along  the  meridian  become 

^  °    -  i 

T,     a.s   =0   ats^s-jands^Sp 

i  =  0 
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1°        -i 
Z     b.s      =0      ats    =    s-|ands    =    Sp 

1=0 


12 


-1 


Z     c.s      =0      ats=Siands=Sp 
i  =  0 


12 

z 

1=0 


2  ,    .    2 


i-2 


Z     c.[(i)(i-l)    +  v(-n^/sin''a   +    i)]s'""      =    0 


at    s    =    s 


and    s    =    s, 
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APPENDIX  II 
CONVERGENCE 

In  the  use  of  the  subdomain  method,  convergence 
in  the  mean  is  desired.   The  main  influence  on  convergence 
is  the  choice  of  trial  functions.   For  assured  convergence, 
the  trial  functions  must  be  complete  and  linearly  indepen- 
dent [49].   The  completeness  property  of  a  set  of  functions 
insures  that  we  can  represent  the  exact  solution  provided 
enough  terms  are  used.   The  power  series  trial  functions 
used  in  this  study  are  complete,  for  example,  so  that  any 
continuous  function  can  be  expanded  in  terms  of  them. 

To  demonstrate  the  rate  of  convergence  for  the 
problem  under  study,  a  series  of  computer  runs  were  con- 
ducted to  determine  the  dynamic  perturbation  response  for 
different  numbers  of  subdomains.   Dynamic  stability  response 
for  eight,  ten,  and  sixteen  subdomains  was  investigated. 
The  results  are  shown  in  Figure  11.   The  ordinate  designa- 
tion implies  percent  differences  with  respect  to  the  best 
approximation. 
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Number  of  Subdomains 


Figure  11.   Convergence  Properties  of  Solution 


APPENDIX  III 


APPENDIX  III 
FLOW  DIAGRAM  OF  THE  COMPUTER  PROGRAM 

For  a  better  understanding  of  the  computer  program, 
a  key  which  may  be  helpful  in  going  from  the  theory  to  the 
program  is  given  below: 


FORTRAN  Name    Theory 


N 

m 

NW 

n 

Q 

Q 

AP 

a 

RO 

P 

E 

E 

V 

V 

F 

^2 

BTT 

B 

STT 

[k] 

DND 

[m] 

PV 

{p} 

PLUG 

{NL} 

X 

T 

Description 

Number  of  subdomains 

Circumferential  wave  number 

Total  number  of  n 

Semi  vertex  angle  of  cone 

Material  density 

Young ' s  modul us 

Poisson's  ratio 

Conical  slant  height 

Boundary  condition  matrix 

Linear  stiffness  matrix 

Mass  matrix 

Load  vector 

Nonlinear  stiffness  vectors 

Nondimensional  time 
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FORTRAN  Name    Theory 


U,V,Y 

n      n      n 

P 

P 

PSD 

Po 

PRMT(2) 

"^2 

PRMT(4) 

e 

ORAN(I) 

-- 

DHPCG 

_  _ 

DGELG 


Descriptions 

Nondimensi onal  displacements 

Axi  al  load 

Pressure    level    of  Heaviside  pulse 

Upper    time    limit 

Error    bound 

Online  core  storage  matrices 

Hamming's  predictor-corrector 
integration  subroutine 

Gauss-elimination  subroutine 


A  flow  diagram  of  the  computer  program  is  illustrated 
in  Figure  12. 
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APPENDIX  IV 
FORTRAN  IV  SOURCE  PROGRAM 

The  computer  program  developed  to  provide  the 
numerical  solution  to  this  dissertation  problem  is  listed 
below.   The  program  requires  182,000  bytes  computer  core 
storage  space.   A  typical  nonlinear  dynamic  stability  run 
executes  in  approximately  in  seven  minutes.   The  compile 
time  is  three  seconds. 
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